#!/bin/bash
set -e

while getopts  "i:p:" opts
do
        case  $opts  in
        
        i)
                        interval=$OPTARG
                        
                ;;
		p)
                        out_prefix=$OPTARG
                        
                ;;
		\?)

				echo `basenaem $0` [-i interval] [-p out_prefix]
				exit 1
				;;
        esac
done
shift $(($OPTIND - 1))


#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------


echo out prefix is: $out_prefix


###--------------- argument may be changed ---------------###


# genome_name=b37.fa
echo using $genome_name as the reference genome
# genome_assembly=b37
echo genome assembly is: $genome_assembly
# dbsnp_version=138
echo dbsnp version is: $dbsnp_version

# data_thread_num=8
# cpu_thread_num=4

# java_memory=16g
echo java memory: $java_memory

# snpeff_db_version=GRCh37.75
echo snpeff database: $snpeff_db_version

# header=\
# @HD\tVN:1.4\tGO:none\tSO:coordinate

# read_group=\
# @RG\\tID:${sample_name}\\tPL:ILLUMINA\\tSM:$sample_name
# For gatk, Each read group must contain the platform (PL) and sample (SM) tags.
# For the platform value, we currently support 454, LS454, Illumina, Solid, ABI_Solid, and CG (all case-insensitive).
# Each read in the file must be associated with exactly one read group.


vcf_path=${data_path}/vcf/gatk/

# there is no space between "=" and "\"
# ref_genome=\
# ${data_path}/ref/b37/$genome_name

echo ref genome is: $ref_genome


reads_seq_1=$1
reads_seq_2=$2

echo reads are: $reads_seq_1 and $reads_seq_2


###--------------- end argument may be changed ---------------###


if [ -n "$interval" ]; then

echo
echo
echo gatk RealignerTargetCreator
java -Xmx$java_memory -jar $gatk \
	-T RealignerTargetCreator \
	-R $ref_genome \
	-I $out_prefix.sort.bam \
	-o $out_prefix.realn.intervals \
	-log $out_prefix.realn.log \
	-nt $data_thread_num \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	-L $interval
	
else

echo
echo
echo gatk RealignerTargetCreator
java -Xmx$java_memory -jar $gatk \
	-T RealignerTargetCreator \
	-R $ref_genome \
	-I $out_prefix.sort.bam \
	-o $out_prefix.realn.intervals \
	-log $out_prefix.realn.log \
	-nt $data_thread_num \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf
	
fi
	
echo
echo
# indel realign
echo gatk IndelRealigner
java -Xmx$java_memory -jar $gatk \
	-T IndelRealigner \
	-R $ref_genome \
	-targetIntervals $out_prefix.realn.intervals \
	-I $out_prefix.sort.bam \
	-o $out_prefix.realn.bam \
	-log $out_prefix.realn2.log \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf